function [gbest,gbestval,fitcount,fitness,gbestval_history]= FCO_func(fhd,Dimension,PopulationSize,MaxIter,lowerbound,upperbound,pos,bool_Plot,varargin)

rand('state',sum(100*clock));
mi = MaxIter;
ps = PopulationSize;
D = Dimension;
[pub, plb, vub, vlb] = SetBound(D,ps,lowerbound,upperbound);

%  Parameters of the Algorithm //////////////////////////////////////
cc=[2 2];  % For PSO
iwt=0.9-(1:mi).*(0.5./mi); % For PSO
G0 = 0.28; % For SGSA
epsilon = 1e-3; % For SGSA
g = 0.13; % For QGSA
NumSC = 4; % For Social Class
%  Parameters of the Algorithm //////////////////////////////////////

%  Settings of the Algorithm //////////////////////////////////////
% If ps = 50, NumSC = 4
% size_each_level = 12    12    12    14
% lb_level = 1    13    25    37
% ub_level = 12    24    36    50
size_each_level = repmat(floor(ps/NumSC),1,NumSC);
size_each_level(NumSC) = ps - (size_each_level(1)*(NumSC-1));
lb_level(1,1) = 1;
ub_level(1,1) = size_each_level(1);
for l = 2 : NumSC
    lb_level(l,1) = sum(size_each_level(1,1:l-1))+1;
    ub_level(l,1) = sum(size_each_level(1,1:l));
end
%  Settings of the Algorithm //////////////////////////////////////

vel = zeros(ps,D); % initialize the velocity of the particles
acc = zeros(ps,D); % initialize the acceleration of the particles

fitness = feval(fhd,pos', varargin{:});
fitcount = ps;
pbest = pos;
pbestval = fitness;
[gbest_fit, gbest_id] = min(pbestval);
gbest = pbest(gbest_id, :);
gbestrep = repmat(gbest, ps, 1);
cbest_fit = min(fitness); % Current best
cworst_fit = max(fitness); % Current worst
gbestval = gbest_fit;
gbestval_history = zeros(fix(mi/500),2);

if(bool_Plot)
    h = figure(8);
    haxes = plot( 0 , 0 );
    XArray = [];
    YArray = [];
    title('FCO');
end
%     h = figure(8);
%     haxes = plot( 0 , 0 );
%     XArray = [];
%     YArray = [];
%     title('FCO');
dcc1 = ones(ps,D).*2;
dcc2 = ones(ps,D).*2;
diwt = ones(ps,D).*0.5;
leader = zeros(NumSC,D);
leader_id = zeros(NumSC,1);
temp_pos = zeros(ps,D);

for iter = 2 : mi    
%     a=2-iter*((2)/mi);

%     if(mod(iter,5)==0)
%         XArray = [ XArray iter]; 
%         YArray = [ YArray sum(sort_mass(ps-fix(ps*0.1):ps,1))];
%         set( haxes , 'XData' , XArray , 'YData' , YArray );
%         drawnow
%     end
    
    rank = [sort_mass(:,2) (1:ps)'];
    rank = sortrows(rank,1);
    rank = rank(:,2);
    leader = zeros(NumSC,D);
    for l = 1 : NumSC
        leader(l,:) = pos(sort_mass(ub_level(l),2),:);
    end
    leaderrep = repmat(leader(1,:),size_each_level(1),1);
    for l = 2 : NumSC
        leaderrep = [leaderrep; repmat(leader(l,:),size_each_level(l),1)];
    end
    
    if(mod(iter,30) == 0)
        dcc1 = zeros(ps,1);
        dcc2 = zeros(ps,1);
        diwt = zeros(ps,1);
        for n = 1 : ps
            if(rank(n) <= ub_level(1))
%                 dcc(n,:) = repmat(16*((mi-iter)/mi)+2,1,D);
                dcc1(n,1) = 2+((iter/mi)^4)*16;
                dcc2(n,1) = 2+(1-((iter/mi)^4))*16;
                diwt(n,1) = 0.1*(1-((iter/mi)^4));
            elseif(rank(n) <= ub_level(2))
%                 dcc(n,:) = repmat(8*((mi-iter)/mi)+2,1,D);
                dcc1(n,1) = 2+((iter/mi)^4)*8;
                dcc2(n,1) = 2+(1-((iter/mi)^4))*8;
                diwt(n,1) = 0.2*(1-((iter/mi)^4));
            elseif(rank(n) <= ub_level(3))
%                 dcc(n,:) = repmat(4*((mi-iter)/mi)+2,1,D);
                dcc1(n,1) = 2+((iter/mi)^4)*4;
                dcc2(n,1) = 2+(1-((iter/mi)^4))*4;
                diwt(n,1) = 0.3*(1-((iter/mi)^4));
            else
                dcc1(n,1) = 2;
                dcc2(n,1) = 2;
                diwt(n,1) = 0.4;
            end
        end
        dcc1 = repmat(dcc1,1,D);
        dcc2 = repmat(dcc2,1,D);
        diwt = repmat(diwt,1,D);
    end
%     acc = dcc1.*rand(ps,D).*(leaderrep-pos) + dcc2.*rand(ps,D).*(gbestrep-pos);
    acc = cc(1).*rand(ps,D).*(pbest-pos) + cc(2).*rand(ps,D).*(leaderrep-pos);
%     if(iter > fix(mi/2))
% %         acc = cc(1).*rand(ps,D).*(pbest-pos) + cc(2).*rand(ps,D).*(gbestrep-pos);
%         acc = cc(1).*rand(ps,D).*(pbest-pos) + cc(2).*rand(ps,D).*(leaderrep-pos);
%     else
%         acc = cc(1).*rand(ps,D).*(leaderrep-pos) + cc(2).*rand(ps,D).*(gbestrep-pos);
%     end
%     acc = dcc1.*rand(ps,D).*(pbest-pos) + dcc2.*rand(ps,D).*(gbestrep-pos);
%     vel = iwt(iter).*vel + acc;
%     vel = diwt.*vel + acc;
    vel = 0.3.*vel + acc;
%     vel = acc;
        
    vel = ((vel>=vlb)&(vel<=vub)).*vel+(vel<vlb).*vlb+(vel>vub).*vub;
    pos = pos+vel;
    pos=((pos>=plb)&(pos<=pub)).*pos...
        +(pos<plb).*(plb+0.25.*(pub-plb).*rand(ps,D))+(pos>pub).*(pub-0.25.*(pub-plb).*rand(ps,D));
    
    fitness = feval(fhd,pos', varargin{:});
    fitcount = fitcount+ps;
    tmp = (pbestval < fitness);
    temp=repmat(tmp',1,D);
    pbest=temp.*pbest+(1-temp).*pos;
    pbestval=tmp.*pbestval+(1-tmp).*fitness; % update the pbest    
    [gbest_fit, gbest_id] = min(pbestval);
    gbest = pbest(gbest_id, :);
    gbestrep = repmat(gbest, ps, 1);
    cbest_fit = min(fitness); % Current best
    cworst_fit = max(fitness); % Current worst
    gbestval = gbest_fit;
        
    if(bool_Plot & mod(iter,10)== 0)
        XArray = [ XArray iter]; 
        YArray = [ YArray gbestval];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
    
    if(mod(iter,500) == 0)
        gbestval_history(fix(iter/500),1) = gbestval;
        gbestval_history(fix(iter/500),2) = iter;
    end
end

end


